library(Seurat)
library(tidyseurat)
library(tidySummarizedExperiment)
library(harmony)
library(SingleR)
library(ggpubr)
library(ggrepel)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

# import ----
c31a1.wd <- 'mission/slc31a1/'

mtx.path <- list.files(c31a1.wd, 'ctrl|porosis',
                       full.names = T, include.dirs = T)

sobj <- mtx.path |>
  str_extract('\\w+$') |>
  set_names(mtx.path, nm = _) |>
  Read10X(strip.suffix = T) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

sobj |>
  VlnPlot('mito.ratio', pt.size = 0) +
  geom_hline(aes(yintercept = 10))

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'ctrl|osteo.+'))

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident','group'), regress.cc = T)

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca_main')

sobj |>
  DimPlot(group.by = 'hpca_main', cols = DiscretePalette(20),
          label = T, repel = T)

sobj |>
  FeaturePlot(c('CSF1R','SLC31A1'))

# rds checkpoint -------
sobj |> write_rds('mission/slc31a1/dong24nrf2.rds')

sobj.mmt <- sobj |>
  filter(hpca_main %in% c('Macrophage', 'Monocyte', 'Tissue_stem_cells'))

sobj.mmt <- sobj.mmt |>
  quick_process_seurat(batch = c('orig.ident','group'), regress.cc = T,
                       skip_norm = T)

sobj.mmt |>
  DimPlot(group.by = 'hpca_main')

sobj.mmt <- sobj.mmt |>
  mark_cell_type_singler(hpca, fine_label = T, new_label = 'hpca_fine')

sobj.mmt |>
  DimPlot(group.by = 'hpca_fine', cols = 'Paired')

sobj.mmt |>
  VlnPlot('CSF1R')

# CSF1R+ Macro/mono as pre-osteoclast --------
sobj.pre.osc <- sobj.mmt |>
  filter(seurat_clusters %in% c(4,5))

sobj.pre.osc |>
  filter(str_detect(hpca_main, 'Macro')) |>
  FeaturePlot('SLC31A1', split.by = 'group')

key.pre.osc <- sobj.pre.osc |>
  get_abundance_sc_long('SLC31A1') |>
  mutate(group = str_extract(.cell, 'ctrl|osteoporosis'))

key.pre.osc |>
  summarise(mean_expr = mean(expm1(.abundance_RNA)),
            pos.ratio = sum(.abundance_RNA > 0)/ n(),
          .by = group)

sobj.pre.osc |>
  FeaturePlot('SLC31A1', split.by = 'group')

sobj.pre.osc |>
  VlnPlot('SLC31A1', group.by = 'group', pt.size = 0) +
  xlab('Group') + NoLegend()

key.pre.osc |>
  ggplot(aes(group, .abundance_RNA, fill=group)) +
  geom_violin(scale = 'width') +
  theme_pubr()

sobj.pre.osc |>
  FindMarkers(features = 'SLC31A1',
              group.by = 'group', ident.1 = 'osteoporosis')

# BMSC -----
sobj.bmsc <- sobj.mmt |>
  filter(str_starts(hpca_fine, 'Tissue'))

sobj.bmsc |>
  FeaturePlot('SLC31A1', split.by = 'group')

key.bmsc <- sobj.bmsc |>
  get_abundance_sc_long('SLC31A1') |>
  mutate(group = str_extract(.cell, 'ctrl|osteoporosis'))

key.bmsc |>
  summarise(mean_expr = mean(expm1(.abundance_RNA)),
            pos.ratio = sum(.abundance_RNA > 0)/ n(),
            .by = group)

key.bmsc |>
  ggplot(aes(group, .abundance_RNA, fill=group)) +
  geom_violin(scale = 'width') +
  theme_pubr(legend = 'none') +
  labs(x = 'Group', y = 'Expression Level', title = 'SLA31A1')

sobj.bmsc |>
  VlnPlot('SLC31A1', group.by = 'group', pt.size = 0) +
  xlab('Group') + NoLegend()

sobj.bmsc |>
  FindMarkers(features = 'SLC31A1',
              group.by = 'group', ident.1 = 'osteoporosis')

sobj.3t <- sobj.mmt |>
  filter(str_starts(hpca_fine, 'Tissue') | seurat_clusters %in% c(4,5)) |>
  mutate(manual_fine = case_when(
    str_detect(hpca_fine, 'Mono') ~ 'c-Fms+ Monocyte',
    str_detect(hpca_fine, 'Macro') ~ 'c-Fms+ Macrophage',
    .default = 'BMSC'),
    manual_main = ifelse(manual_fine == 'BMSC', 'BMSC', 'Osteoclast precursor'))

sobj.3t |>
  VlnPlot('SLC31A1', group.by = 'manual_fine', pt.size = 0) +
  NoLegend() + xlab('Cell type')

sobj.3t |>
  VlnPlot('SLC31A1', group.by = 'manual_main', pt.size = 0) +
  NoLegend() + xlab('Cell type')

sobj.3t |>
  FindMarkers(features = 'SLC31A1', group.by = 'manual_main',
              ident.1 = 'BMSC')

sobj.3t |>
  FindMarkers(features = 'SLC31A1', group.by = 'manual_fine',
              ident.1 = 'c-Fms+ Macrophage', ident.2 = 'c-Fms+ Monocyte')
